import FE
from matplotlib import pyplot as plt




def u_true(x):
    return x**4 * (1 - x)**4


def f(x):
    return 24 * (x**4 + 16 * x**3 * (x - 1) + 36 * x**2 *
                 (x - 1)**2 + 16 * x * (x - 1)**3 + (x - 1)**4)

n = 16
msh = FE.Mesh(n, interval=[0, 1])
fe = FE.FE(msh)
fe.solve(f, left_value=0, right_value=0)
u_values = fe.u[::fe.n_dofs_per_node]  # 因为每个节点有两个自由度，这里提取函数值
l2_norm = fe.compute_l2_norm(u)
grad_norm = fe.compute_gradient_l2_norm(u)
hessian_norm = fe.compute_hessian_l2_norm(u)
error_l2 = fe.compute_l2_error(u, true_u)

print(f"L2范数: {l2_norm:.6e}")
print(f"一阶导数L2范数: {grad_norm:.6e}")
print(f"二阶导数L2范数: {hessian_norm:.6e}")
print(f"L2误差: {error_l2:.6e}")

# 绘制解的图像
plt.plot(msh.mesh, u_values, marker='o', label='FEM')
plt.plot(msh.mesh, u_true(msh.mesh), marker='*', label="Exact")
plt.xlabel('x')
plt.ylabel('u(x)')
plt.legend()
plt.title('Solution of the Finite Element Method')
plt.grid(True)
plt.show()
